Species change at pixel level

library(tidyverse)
library(terra)
library(sf)
library(rnaturalearth)
library(here)
library(iucnredlist) ### remotes::install_github("IUCN-UK/iucnredlist")
library(janitor)
library(cowplot)

This document investigates the Sorensen Similarity Index for the bird species in the Lawler paper at both the species and trait level. Visualizations are at the pixel level. Categorical trait data is from AVONET.

# Functions 
source(here("common_fxns.R"))

# Read in Lawler data
spp_mapfiles <- list.files(here_anx('/data_lawler_2020/species_projections'), 
                           full.names = TRUE, recursive = TRUE, pattern = '.tif')
spp_all_df <- data.frame(f = spp_mapfiles,
                         spp = basename(spp_mapfiles) %>%
                           str_remove('(_gf85|_in85|_mc85|_bias).+') %>%
                           str_replace_all('_', ' '),
                         tax = basename(dirname(spp_mapfiles)))

# Read in AVONET dataset (reading in BirdLife tab, which should match with IUCN names)
avonet <- read_csv(here("data-laa", "AVONET-birdlife-tab-supplementary-dataset-1.csv")) %>% 
  clean_names()

# Filter for just Lawler birds
birds_all_df <- spp_all_df %>% 
  filter(tax == "Birds")

birds_all <- birds_all_df %>% 
  group_by(spp, tax) %>% 
  summarize(files = n())

Gather maps and aggregate

Following the method in 2_combine_spp_scenarios.qmd: For all the Lawler birds species (301) gather all future scenario maps and aggregate, summing and dividing by 3 to make a “probability” map for future distributions, and crop to just the lower 48 US states (and a little Canada and Mexico). Then crop the historical maps to the same extent.

Aggregate future projections

For each species in the sample, aggregate the three future scenarios into a single map. Crop to the lower continental US states.

### set up extent for continental US
us_extent <- ext(c(xmin = -130, xmax = -65, ymin = 23, ymax = 50))

### Aggregate over the vector of species names
spp_vec <- birds_all_df$spp %>% unique()

future_rast <- lapply(spp_vec, FUN = agg_scenario, df = birds_all_df) %>%
  setNames(spp_vec) %>%
  rast()

x <- future_rast[[1]]
#plot(x, main = names(future_rast)[1])

crop historical maps

hist_rast <- lapply(spp_vec, FUN = agg_scenario, df = birds_all_df, 
                    scenario = 'historical') %>%
  setNames(spp_vec) %>%
  rast()

x <- hist_rast[[1]]
#plot(x, main = names(hist_rast)[1])

Generate jurisdiction map

Let’s generate a simple jurisdiction map - US states, and Canada and Mexico.

states_sf <- ne_states(returnclass = 'sf', country = c('Mexico', 'Canada', 'United States of America')) 

state_id_df <- states_sf %>%
  st_drop_geometry() %>%
  mutate(juris = ifelse(adm0_a3 == 'USA', name, admin)) %>%
  dplyr::select(id = gn_id, juris)
### plot(states_sf %>% select(iso_a2))
states_rast <- rasterize(states_sf, hist_rast, field = 'gn_id')
#plot(states_rast)

# Function to get pixels instead of jurisdictions
extract_jurisdictions_pixels <- function(spp, future_r, hist_r, juris_r) {
  ### spp <- spp_vec[1]
  ### future_r = future_rast; hist_r = hist_rast; juris_r = states_rast
  s_f_r <- future_r[[names(future_r) == spp]]
  s_h_r <- hist_r[[names(hist_r) == spp]]
  s_df <- data.frame(values(s_h_r),
                     values(s_f_r),
                     values(juris_r)) %>%
    setNames(c('historical', 'future', 'id')) %>%
    mutate(pixel_id = row_number()) %>% 
    filter(!is.na(id)) %>%
    filter(!is.na(historical) | !is.na(future)) #%>%
    #group_by(id) %>%
    #summarize(n_historical = sum(historical > .25),
              #n_future = sum(future > .5))

  return(s_df)
}

Rough pass at classification?

For each spp let’s create a dataframe of jurisdictions and cell values…

state_spp_list <- lapply(spp_vec, FUN = extract_jurisdictions_pixels,
                         future_r = future_rast, 
                         hist_r = hist_rast,
                         juris_r = states_rast) %>%
  setNames(spp_vec) %>%
  bind_rows(.id = 'spp') %>%
  left_join(state_id_df, by = 'id') #%>%
  ### note, Canada and Mexico have multiple jurisdictions (provinces and states)
  ### so sum up all those instances
  #group_by(spp, juris) %>%
  #summarize(n_historical = sum(n_historical),
            #n_future = sum(n_future), 
            #.groups = 'drop')

Pixel level analysis

Now that we have each species at the pixel level (pixel id), how many species are in each pixel for historical and future?

pixel_spp <- state_spp_list %>% 
  group_by(pixel_id, juris) %>% 
  summarize(n_spp_historical = sum(historical > .25),
            n_spp_future = sum(future > .5)) %>% 
  mutate(n_spp_diff = n_spp_future - n_spp_historical)

# Vast majority of pixels are losing total # of species (simply from numbers perspective)
ggplot(data = pixel_spp) +
  geom_histogram(aes(x = n_spp_diff),
                 bins = 50) +
  labs(title = "Bird species gain from historic to future by pixel",
       x = "species change")

Sorensen Similarity Index

# Test out with a single pixel (141325)
pxl_141325 <- state_spp_list %>% 
  filter(pixel_id == 141325) %>% 
  # Assign future pixels > 0.5 to spp prescence 
  mutate(future = case_when(future > 0.5 ~ 1,
                           TRUE ~ 0))

pxl_si <- pxl_141325 %>% 
  group_by(pixel_id) %>% 
  summarize(shared = sum(historical == 1 & future == 1),
            historical_spp = sum(historical),
            future_spp = sum(future),
            ssi = (2 *shared)/(historical_spp + future_spp))

# Try all pixels
pxl_ssi <- state_spp_list %>% 
  # Assign future pixels > 0.5 to spp prescence 
  mutate(future = case_when(future > 0.5 ~ 1,
                           TRUE ~ 0)) %>% 
  group_by(pixel_id, juris) %>% 
  summarize(shared = sum(historical == 1 & future == 1),
            historical_spp = sum(historical),
            future_spp = sum(future),
            ssi = (2 *shared)/(historical_spp + future_spp))

# Distribution of Sorensen Similarity Index for bird species
# Mean is 0.40 (moderate similarity between historic and future)
ggplot(data = pxl_ssi) +
  geom_histogram(aes(x = ssi),
                 bins = 50) +
  geom_vline(xintercept = mean(pxl_ssi$ssi)) +
  annotate("text", label = "mean = 0.40", x = 0.30, y = 15000) +
  #facet_wrap(~juris) +
  labs(title = "Sorensen Similarity Index btwn historic & future by pixel (n = 152,594)",
       x = "SSI")

# SSI stats by jurisdiction
ssi_juris <- pxl_ssi %>% 
  group_by(juris) %>% 
  summarize(mean_ssi = mean(ssi),
            median_ssi = median(ssi),
            max_ssi = max(ssi),
            min_ssi = min(ssi)) %>% 
  ungroup()

ssi_juris %>% 
  mutate(juris = fct_reorder(juris, mean_ssi, .desc = TRUE)) %>%
  ggplot() +
  geom_col(aes(x = juris, y = mean_ssi)) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  labs(y = "mean pixel SSI",
       x = "jurisdiction",
       title = "Bird species mean pixel SSI by jurisdiction")

Map the SSI by pixel

# Create a blank raster with same extent as x = hist_rast[[1]] 
id_rast <- rast(x)

# Assign IDs and then convert to sf
values(id_rast ) <- 1:ncell(x)

id_rast_sf <- id_rast %>% 
  stars::st_as_stars() %>% 
  st_as_sf() %>% 
  rename("pixel_id" = "Acrocephalus familiaris")

# Merge in the SSI
id_rast_sf_ssi <- inner_join(id_rast_sf, pxl_ssi)

ggplot() +
  geom_sf(data = id_rast_sf_ssi %>% 
            filter(juris == "California"), aes(geometry = geometry, fill = ssi, color = ssi), size = 0) + 
  scale_fill_viridis_c(limits = c(0,1)) + 
  scale_color_viridis_c(limits = c(0,1)) +
  labs(title = "Birds in CA Sorensen Similiarity Index")

ggplot() +
  geom_sf(data = id_rast_sf_ssi, aes(geometry = geometry, color = ssi, fill = ssi)) +
  scale_color_viridis_c() + 
  scale_fill_viridis_c() + 
  labs(title = "Birds Sorensen Similiarity Index")

In general, SSI values are as follows (https://www.statology.org/how-to-calculate-and-interpret-the-sorensen-dice-index/):

• 0.80-1.00: Very high similarity (communities share most species)
• 0.60-0.79: High similarity (substantial species overlap)
• 0.40-0.59: Moderate similarity (some shared species)
• 0.20-0.39: Low similarity (few shared species)
• 0.00-0.19: Very low similarity (minimal species overlap)

Species numbers

historic_spp_num <- ggplot() +
  geom_sf(data = id_rast_sf_ssi, aes(geometry = geometry, fill = historical_spp, color = historical_spp)) + 
  scale_fill_viridis_c(option = "plasma", limits = c(1,180)) + 
  scale_color_viridis_c(option = "plasma", limits = c(1,180)) +
  geom_sf(data = id_rast_sf_ssi %>% filter(historical_spp == 0), aes(geometry = geometry), fill = "gray40", color = "grey40") + 
  labs(title = "Historic number of bird species",
       subtitle = "(0 species in gray)",
       fill = "count",
       color = "count")

future_spp_num <- ggplot() +
  geom_sf(data = id_rast_sf_ssi, aes(geometry = geometry, fill = future_spp, color = future_spp)) + 
  scale_fill_viridis_c(option = "plasma", limits = c(1,180)) + 
  scale_color_viridis_c(option = "plasma", limits = c(1,180)) +
  geom_sf(data = id_rast_sf_ssi %>% filter(future_spp == 0), aes(geometry = geometry), fill = "gray40", color = "grey40") + 
  labs(title = "Future number of bird species",
       subtitle = "(0 species in gray)",
       fill = "count",
       color = "count")

plot_grid(historic_spp_num + theme(legend.position = 'bottom'), 
          future_spp_num + theme(legend.position = 'bottom'))

Add in Avonet data and investigate trait SSI at pixel level

Let’s start with the birds that easily merge with the Avonet dataset (270)

# See how well the Avonet trait data merges in 
# 270 of 301 have matches - could work with these species for now since they should be able to match with IUCN
birds_avonet <- left_join(birds_all, (avonet %>% rename(spp = species1)))

lawler_avvonet_birds <- birds_avonet %>% 
  drop_na(sequence)

lawler_avvonet_birds_vec <- lawler_avvonet_birds$spp

Now, will add in the bird species that did not merge as easily for a total of 299 (2 of the 301 birds are extinct)

# Read in crosswalk of lawler to avonet species (created with input from gemini)
lawler_avonet_xwlk <- read_csv(here("data-laa", "lawler_avonet_xwlk.csv")) %>% 
  clean_names() %>% 
  dplyr::select(lawler_spp:name_status) %>% 
  drop_na()

# Match the lawler avonet xlk with trait data
# The two extinct species drop out
lawler_avonet_traits_missing <- inner_join(lawler_avonet_xwlk, 
                                   avonet %>% rename(avonet_spp = species1))

# Bind the easily matched birds and the birds that were missing in the merge
lawler_avonet_all <- rbind(
  lawler_avvonet_birds %>% dplyr::select(-tax, -files),
  lawler_avonet_traits_missing %>% dplyr::select(-avonet_spp, -name_status)
)

pxl_spp_list_trait <- inner_join(state_spp_list, lawler_avonet_all) %>% 
    filter(!is.na(id)) %>%
    filter(!is.na(historical) | !is.na(future))

Create trait SSI map

us_spp_trait <- pxl_spp_list_trait %>% 
  dplyr::select(spp:juris, habitat:primary_lifestyle) %>% 
  filter(historical > 0.25 | future > 0.5) %>% 
  mutate(future = case_when(future > 0.5 ~ 1,
                           TRUE ~ 0)) %>% 
  mutate(habitat_density = case_when(habitat_density == 1 ~ "Dense",
                                     habitat_density == 2 ~ "Semi-open",
                                     habitat_density == 3 ~ "Open")) %>% 
  mutate(migration = case_when(migration == 1 ~ "Sedentary",
                               migration == 2 ~ "Partial",
                               migration == 3 ~ "Migratory"))

# Combine the categorical columns into one group column to create unique groups
# Count up these trait groups by pixel
us_spp_trait_groups <- us_spp_trait %>% 
  unite(trait_group, migration, habitat_density, habitat, trophic_level, trophic_niche, primary_lifestyle) %>% 
  group_by(pixel_id, juris, trait_group) %>% 
  summarize(historical_trait_group_count = sum(historical),
            future_trait_group_count = sum(future)) %>% 
  ungroup()

# If we apply the Sorensen Similarity Index to traits, then we care about the presence / absence of a trait group and not the # of species
us_spp_trait_groups_ssi <- us_spp_trait_groups %>% 
  mutate(historical_trait_group_presence = case_when(historical_trait_group_count == 0 ~ 0,
                                                     historical_trait_group_count != 0 ~ 1),
         future_trait_group_presence = case_when(future_trait_group_count == 0 ~ 0,
                                                     future_trait_group_count != 0 ~ 1)) %>% 
  group_by(pixel_id, juris) %>% 
  summarize(shared = sum(historical_trait_group_presence == 1 & future_trait_group_presence == 1),
            historical_trait_groups= sum(historical_trait_group_presence),
            future_trait_groups = sum(future_trait_group_presence),
            trait_group_ssi = (2 *shared)/(historical_trait_groups + future_trait_groups)) %>% 
  ungroup()

# Merge in pixel geography for map
id_rast_sf_us_ssi <- inner_join(id_rast_sf, us_spp_trait_groups_ssi)

# Create map (299 bird species)
ggplot() +
  geom_sf(data = id_rast_sf_us_ssi, aes(geometry = geometry, fill = trait_group_ssi, color = trait_group_ssi), size = 0) + 
  scale_fill_viridis_c(limits = c(0,1)) + 
  scale_color_viridis_c(limits = c(0,1)) +
  labs(title = "Bird trait groups Sorensen Similiarity Index")

What’s the change between trait SSI and species SSI?

If traits are compensating for some of the change in species similarity, then the difference between trait SSI and species SSI should be positive (more similarity between traits than species).

pxl_ssi_dif <- left_join(us_spp_trait_groups_ssi %>% 
                           dplyr::select(pixel_id, trait_group_ssi),
                         pxl_ssi %>% 
                           dplyr::select(pixel_id, ssi)) %>% 
  mutate(ssi_diff = trait_group_ssi - ssi)

# Merge in the pixel geography
id_rast_sf_us_ssi_diff <- inner_join(id_rast_sf, pxl_ssi_dif)

ggplot() +
  geom_sf(data = id_rast_sf_us_ssi_diff, aes(geometry = geometry, fill = ssi_diff, color = ssi_diff), size = 0) + 
  scale_color_distiller(palette = "BrBG", direction = 1, limits = c(-0.4, 0.4)) + 
  scale_fill_distiller(palette = "BrBG", direction = 1, limits = c(-0.4, 0.4)) +
  # Show areas 0 species and trait overlap in pink
  geom_sf(data = id_rast_sf_us_ssi_diff %>% 
            filter(ssi == 0 & trait_group_ssi == 0), 
          aes(geometry = geometry), fill = "pink", color = "pink") +
  # Show areas with full species and trait overlap in navy
  geom_sf(data = id_rast_sf_us_ssi_diff %>% 
            filter(ssi == 1 & trait_group_ssi == 1), 
          aes(geometry = geometry), fill = "navy", color = "navy") +
  labs(title = "Trait SSI - Species SSI",
       subtitle = "(pink is 0 trait and 0 species SSI; navy is 1 trait and 1 species SSI)")

Areas in teal have greater trait similarity than species similarity. Areas in brown have less trait similarity than species similarity.

Let’slook closer at the pixels where the difference is 0

Zeros may be areas where there are 0 species SSI and 0 trait SSI (no similarity btwn historic and future), 1 species SSI and 1 trait SSI (complete similarity btwn historic and future), or nonzero values.

pxl_ssi_diff_0 <- pxl_ssi_dif %>% 
  filter(ssi_diff == 0)

pxl_ssi_diff_0_summary <- pxl_ssi_diff_0 %>% 
  mutate(diff_type = case_when(trait_group_ssi == 0 & ssi == 0 ~ "zero",
                               trait_group_ssi == 1 & ssi == 1 ~ "one",
                               TRUE ~ "non-zero")) %>% 
  group_by(diff_type) %>% 
  summarize(count = n()) %>% 
  ungroup()

zero_trait_species_ssi <- pxl_ssi_diff_0 %>% 
  filter(trait_group_ssi == 0 & ssi == 0 )

zero_trait_species_ssi_pxls <- zero_trait_species_ssi$pixel_id

# Look at the zero ssi pixels

zero_species_ssi <-state_spp_list %>% 
  # Assign future pixels > 0.5 to spp prescence 
 # filter(spp %in% lawler_avvonet_birds_vec) %>% 
  filter(pixel_id %in% zero_trait_species_ssi_pxls) %>% 
  mutate(future = case_when(future > 0.5 ~ 1,
                           TRUE ~ 0)) %>% 
  # Assigning presence / absence for these birds
  filter(historical == 1 | future == 1)

Plotting a few of the North Dakota birds

North Dakota looks to be an area with high species loss. These maps take a closer look at a few of the species there.

Xenus cinereus - Terek sandpiper

This Eurasian species has only been observed 5 times in the lower 48 states. Surprising to see the historical model project such a large range. This is one of the few species moving into North Dakota (although seems unlikely in reality).

plot(hist_rast$`Xenus cinereus`, main = "Xenus cinereus - historical")

plot(future_rast$`Xenus cinereus`, main = "Xenus cinereus - future")

Parus hudsonicus - boreal chickadee

plot(hist_rast$`Parus hudsonicus`, main = "Parus hudsonicus - historical")

plot(future_rast$`Parus hudsonicus`, main = "Parus hudsonicus - future")

Anser rossii - Ross’s Goose

Large projected reduction in range

plot(hist_rast$`Anser rossii`, main = "Anser rossii - historical")

plot(future_rast$`Anser rossii`, main = "Anser rossii - future")

Tympanuchus pallidicinctus - lesser prairie-chicken

plot(hist_rast$`Tympanuchus pallidicinctus`, main = "Tympanuchus pallidicinctus - historical")

plot(future_rast$`Tympanuchus pallidicinctus`, main = "Tympanuchus pallidicinctus - future")

Two extinct species

These both have about 15 historic GBIF records in the 1900s in Hawaii (around when they went extinct).

plot(hist_rast$`Moho nobilis`, main = "Moho nobilis - historical")

plot(future_rast$`Moho nobilis`, main = "Moho nobilis - future")

plot(hist_rast$`Hemignathus obscurus`, main = "Hemignathus obscurus - historical")

plot(future_rast$`Hemignathus obscurus`, main = "Hemignathus obscurus - future")